import numpy as np
import pandas as pd
from scipy import stats
import matplotlib.pyplot as plt
import plotly.express as px
import seaborn as sns
from patsy import dmatrices
from statsmodels.stats.outliers_influence import variance_inflation_factor
from scipy.stats import chi2_contingency
from statsmodels.stats.contingency_tables import mcnemar
from cmh import CMH
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_absolute_error, mean_squared_error, median_absolute_error, r2_score
from sklearn.preprocessing import PolynomialFeatures
from outliers import smirnov_grubbs as grubbs
from neulab.OutlierDetection import DixonTest
from scipy.stats import kstest, shapiro, anderson, cramervonmises, norm
from statsmodels.stats.diagnostic import lilliefors
from numpy.random import randn
from scipy.stats import ttest_ind
import statsmodels.stats.power as smp
%matplotlib inline
%config InlineBackend.figure_format = 'svg'
import warnings
warnings.filterwarnings('ignore')
Данные: Heart Attack dataset - 13 медицинских признаков и колонка о наличии сердечного приступа
Источник: https://www.kaggle.com/datasets/pritsheta/heart-attack
df = pd.read_csv('Heart_Attack.csv')
df.head()
| age | sex | cp | trestbps | chol | fbs | restecg | thalach | exang | oldpeak | slope | ca | thal | target | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 63 | 1 | 3 | 145 | 233 | 1 | 0 | 150 | 0 | 2.3 | 0 | 0 | 1 | 1 |
| 1 | 37 | 1 | 2 | 130 | 250 | 0 | 1 | 187 | 0 | 3.5 | 0 | 0 | 2 | 1 |
| 2 | 41 | 0 | 1 | 130 | 204 | 0 | 0 | 172 | 0 | 1.4 | 2 | 0 | 2 | 1 |
| 3 | 56 | 1 | 1 | 120 | 236 | 0 | 1 | 178 | 0 | 0.8 | 2 | 0 | 2 | 1 |
| 4 | 57 | 0 | 0 | 120 | 354 | 0 | 1 | 163 | 1 | 0.6 | 2 | 0 | 2 | 1 |
df.info()
<class 'pandas.core.frame.DataFrame'> RangeIndex: 303 entries, 0 to 302 Data columns (total 14 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 age 303 non-null int64 1 sex 303 non-null int64 2 cp 303 non-null int64 3 trestbps 303 non-null int64 4 chol 303 non-null int64 5 fbs 303 non-null int64 6 restecg 303 non-null int64 7 thalach 303 non-null int64 8 exang 303 non-null int64 9 oldpeak 303 non-null float64 10 slope 303 non-null int64 11 ca 303 non-null int64 12 thal 303 non-null int64 13 target 303 non-null int64 dtypes: float64(1), int64(13) memory usage: 33.3 KB
for n in ['age', 'trestbps', 'chol', 'thalach', 'oldpeak']:
sns.displot(x=n,
kde=True,
data=df)
sns.displot(data = df,
x='trestbps',
hue='cp',
multiple='fill',
kind='kde')
<seaborn.axisgrid.FacetGrid at 0x7fb4798b0df0>
sns.boxplot(data=df, x='cp', y='trestbps')
<AxesSubplot:xlabel='cp', ylabel='trestbps'>
fig = px.scatter(df, y='cp', x='trestbps')
fig.update_yaxes(type='category', tickvals=df['cp'].tolist())
fig.show()
plt.rcParams['figure.figsize'] = (30, 3)
sns.stripplot(data=df,
y='cp',
x='trestbps',
jitter=False,
s=15, marker="D", linewidth=1, alpha=0.1)
plt.yticks([0, 1, 2, 3])
([<matplotlib.axis.YTick at 0x7fb47907b8b0>, <matplotlib.axis.YTick at 0x7fb47907b130>, <matplotlib.axis.YTick at 0x7fb477e98940>, <matplotlib.axis.YTick at 0x7fb477e98970>], [Text(0, 0, ''), Text(0, 0, ''), Text(0, 0, ''), Text(0, 0, '')])
chol_ = sorted(list(df.chol))
chol = chol_[140:160]
chol.append(chol_[0])
chol.append(chol_[-1])
plt.rcParams['figure.figsize'] = (20, 3)
sns.boxplot(data=chol, orient = 'h')
<AxesSubplot:>
gtest = grubbs.test(chol)
print('Grubbs test', end=' ')
print(*gtest)
print('chol', end=' ')
print(*chol)
print('min', end=' ')
print(*grubbs.min_test_outliers(chol))
print('max', end=' ')
print(*grubbs.max_test_outliers(chol))
Grubbs test 236 236 236 237 239 239 239 239 240 240 240 240 241 242 243 243 243 243 244 244 chol 236 236 236 237 239 239 239 239 240 240 240 240 241 242 243 243 243 243 244 244 126 564 min 126 max 564
sns.boxplot(data=gtest, orient='h')
<AxesSubplot:>
dtest = DixonTest(dataframe=pd.DataFrame(chol), q=95, info=True, autorm=True)
Detected outlier:
0
21 564
print('Dixon test', end=' ')
print(*dtest[0])
Dixon test 236 236 236 237 239 239 239 239 240 240 240 240 241 242 243 243 243 243 244 244 126
sns.boxplot(data=dtest, orient='h')
<AxesSubplot:>
dtest = DixonTest(dataframe=pd.DataFrame(dtest[0]), q=95, info=True, autorm=True)
Detected outlier:
0
20 126
print('Dixon test', end=' ')
print(*dtest[0])
print('chol', end=' ')
print(*chol)
Dixon test 236 236 236 237 239 239 239 239 240 240 240 240 241 242 243 243 243 243 244 244 chol 236 236 236 237 239 239 239 239 240 240 240 240 241 242 243 243 243 243 244 244 126 564
sns.boxplot(data=dtest, orient='h')
<AxesSubplot:>
trestbps = df[['trestbps']].copy()
trestbps.index = range(trestbps.shape[0])
missing_indexes = np.random.choice(trestbps.shape[0], int(trestbps.shape[0] * 0.1))
trestbps.loc[:, 'miss'] = trestbps.loc[:, 'trestbps']
trestbps.loc[missing_indexes, 'miss'] = np.nan
mean = trestbps['miss'].fillna(trestbps['miss'].mean()).copy()
print('RMSE: ' + str(np.sqrt(np.mean((mean[missing_indexes] - trestbps.loc[missing_indexes, 'trestbps'])**2))))
RMSE: 14.606377063125159
median = trestbps['miss'].fillna(trestbps['miss'].median()).copy()
print('RMSE: ' + str(np.sqrt(np.mean((median[missing_indexes] - trestbps.loc[missing_indexes, 'trestbps'])**2))))
RMSE: 14.704874475266129
seed=42
np.random.seed(seed)
#функция для генерации эталонного распределения
def generate_sample(e, sigma):
x = np.linspace(e - 3*sigma, e + 3*sigma, 1000)
x_ = np.linspace(-3, 3, 1000)
y = norm.cdf(x_)
sns.lineplot(x=x, y=y)
N = {(0, 1, 100) : norm.rvs(loc=0, scale=1, size=100, random_state=seed),
(0, 1, 1000) : norm.rvs(loc=0, scale=1, size=1000, random_state=seed),
(7, 3, 100) : norm.rvs(loc=7, scale=3, size=100, random_state=seed),
(7, 3, 1000) : norm.rvs(loc=7, scale=3, size=1000, random_state=seed)}
for e, sigma, n in [(0, 1, 100), (0, 1, 1000), (7, 3, 100), (7, 3, 1000)]:
sns.ecdfplot(N[(e, sigma, n)])
generate_sample(e, sigma)
plt.show()
for e, sigma, n in [(0, 1, 100), (0, 1, 1000), (7, 3, 100), (7, 3, 1000)]:
stats.probplot(N[(e, sigma, n)], dist="norm", plot=plt)
plt.show()
Реализация только на R. 
for e, sigma, n in [(0, 1, 100), (0, 1, 1000), (7, 3, 100), (7, 3, 1000)]:
print(f'N({e}, {sigma}), n = {n}:', kstest(N[(e, sigma, n)], stats.norm.cdf))
N(0, 1), n = 100: KstestResult(statistic=0.10357070563896065, pvalue=0.21805553378516185) N(0, 1), n = 1000: KstestResult(statistic=0.017327787320720822, pvalue=0.9196626608357358) N(7, 3), n = 100: KstestResult(statistic=0.9348327888393462, pvalue=5.05808847359865e-119) N(7, 3), n = 1000: KstestResult(statistic=0.9436852219799222, pvalue=0.0)
for e, sigma, n in [(0, 1, 100), (0, 1, 1000), (7, 3, 100), (7, 3, 1000)]:
print(f'N({e}, {sigma}), n = {n}:', shapiro(N[(e, sigma, n)]))
N(0, 1), n = 100: ShapiroResult(statistic=0.9898834228515625, pvalue=0.6551706790924072) N(0, 1), n = 1000: ShapiroResult(statistic=0.9986080527305603, pvalue=0.6264819502830505) N(7, 3), n = 100: ShapiroResult(statistic=0.9898834228515625, pvalue=0.6551706790924072) N(7, 3), n = 1000: ShapiroResult(statistic=0.9986080527305603, pvalue=0.6264819502830505)
for e, sigma, n in [(0, 1, 100), (0, 1, 1000), (7, 3, 100), (7, 3, 1000)]:
print(f'N({e}, {sigma}), n = {n}:', anderson(N[(e, sigma, n)]))
N(0, 1), n = 100: AndersonResult(statistic=0.25343395875111696, critical_values=array([0.555, 0.632, 0.759, 0.885, 1.053]), significance_level=array([15. , 10. , 5. , 2.5, 1. ])) N(0, 1), n = 1000: AndersonResult(statistic=0.3474697767348971, critical_values=array([0.574, 0.653, 0.784, 0.914, 1.088]), significance_level=array([15. , 10. , 5. , 2.5, 1. ])) N(7, 3), n = 100: AndersonResult(statistic=0.25343395875111696, critical_values=array([0.555, 0.632, 0.759, 0.885, 1.053]), significance_level=array([15. , 10. , 5. , 2.5, 1. ])) N(7, 3), n = 1000: AndersonResult(statistic=0.3474697767348971, critical_values=array([0.574, 0.653, 0.784, 0.914, 1.088]), significance_level=array([15. , 10. , 5. , 2.5, 1. ]))
for e, sigma, n in [(0, 1, 100), (0, 1, 1000), (7, 3, 100), (7, 3, 1000)]:
print(f'N({e}, {sigma}), n = {n}:', cramervonmises(N[(e, sigma, n)], stats.norm.cdf))
N(0, 1), n = 100: CramerVonMisesResult(statistic=0.1549979480315905, pvalue=0.3751758581771183) N(0, 1), n = 1000: CramerVonMisesResult(statistic=0.06313990902074608, pvalue=0.7943350162473933) N(7, 3), n = 100: CramerVonMisesResult(statistic=31.500937472206292, pvalue=1.2769140411705848e-08) N(7, 3), n = 1000: CramerVonMisesResult(statistic=319.7856715981252, pvalue=1.0722391730055847e-07)
for e, sigma, n in [(0, 1, 100), (0, 1, 1000), (7, 3, 100), (7, 3, 1000)]:
print(f'N({e}, {sigma}), n = {n}:', lilliefors(N[(e, sigma, n)]))
N(0, 1), n = 100: (0.051776473605976925, 0.7370142762533108) N(0, 1), n = 1000: (0.021396139137066617, 0.4050632259039718) N(7, 3), n = 100: (0.05177647360597731, 0.7370142762532998) N(7, 3), n = 1000: (0.02139613913706656, 0.40506322590397575)
sns.displot(x='trestbps', kde=True, aspect=2.5, data=df)
<seaborn.axisgrid.FacetGrid at 0x7fb4763203a0>
#функция для генерации эталонного распределения
def generate_sample(e, sigma):
x = np.linspace(e - 3*sigma, e + 3*sigma, 1000)
x_ = np.linspace(-3, 3, 1000)
y = norm.cdf(x_)
sns.lineplot(x=x, y=y)
trestbps_50 = df.trestbps[0:50]
plt.figure(figsize=(18,7))
plt.subplot(1, 3, 1)
sns.ecdfplot(trestbps_50, color = 'orange')
plt.xlabel('ecdf')
plt.subplot(1, 3, 2)
generate_sample(135, 18)
sns.ecdfplot(trestbps_50, color = 'orange')
plt.xlabel('ecdf, N(135, 18)')
plt.subplot(1, 3, 3)
generate_sample(135, 18)
plt.xlabel('N(135, 18)')
plt.show()
stats.probplot(trestbps_50, dist="norm", plot=plt)
plt.show()
Метод огибающих 
print(kstest(trestbps_50, stats.norm.cdf))
print(shapiro(trestbps_50))
print(anderson(trestbps_50))
print(cramervonmises(trestbps_50, stats.norm.cdf))
print('Lilliefors', lilliefors(trestbps_50))
KstestResult(statistic=1.0, pvalue=0.0) ShapiroResult(statistic=0.9752922654151917, pvalue=0.37441587448120117) AndersonResult(statistic=0.5554230489995859, critical_values=array([0.538, 0.613, 0.736, 0.858, 1.021]), significance_level=array([15. , 10. , 5. , 2.5, 1. ])) CramerVonMisesResult(statistic=16.666666666666664, pvalue=5.220809340400479e-10) Lilliefors (0.10271428206231425, 0.21230434604399845)
trestbps = df.trestbps
plt.figure(figsize=(18,7))
plt.subplot(1, 3, 1)
sns.ecdfplot(trestbps, color = 'orange')
plt.xlabel('ecdf')
plt.subplot(1, 3, 2)
generate_sample(135, 18)
sns.ecdfplot(trestbps, color = 'orange')
plt.xlabel('ecdf, N(135, 18)')
plt.subplot(1, 3, 3)
generate_sample(135, 18)
plt.xlabel('N(135, 18)')
plt.show()
stats.probplot(trestbps, dist="norm", plot=plt)
plt.show()
Метод огибающих 
print(kstest(trestbps, stats.norm.cdf))
print(shapiro(trestbps))
print(anderson(trestbps))
print(cramervonmises(trestbps, stats.norm.cdf))
print("Lilliefors", lilliefors(trestbps))
KstestResult(statistic=1.0, pvalue=0.0) ShapiroResult(statistic=0.965917706489563, pvalue=1.4579997014152468e-06) AndersonResult(statistic=2.5716952861943128, critical_values=array([0.569, 0.648, 0.777, 0.906, 1.078]), significance_level=array([15. , 10. , 5. , 2.5, 1. ])) CramerVonMisesResult(statistic=101.0, pvalue=0) Lilliefors (0.10194559488884891, 0.0009999999999998899)
def student_test(alpha, data1, data2):
# Гипотеза о равенстве
stat, p_value = ttest_ind(data1, data2)
if p_value > alpha:
print('Принять гипотезу о равенстве')
else:
print('Отклонить гипотезу о равенстве')
print("Мощность критерия", smp.ttest_power(np.mean(data1)/np.std(data1),
nobs=len(data1), alpha = alpha, alternative='two-sided'))
print()
# Гипотеза о большем значении среднего первой выборке'
stat, p_value = ttest_ind(data1, data2, alternative = 'less')
if p_value > alpha:
print('Принять гипотезу о большем значении среднего первой выборки')
else:
print('Отклонить гипотезу о большем значении среднего первой выборки')
print("Мощность критерия", smp.ttest_power(np.mean(data1)/np.std(data1),
nobs=len(data1), alpha = alpha, alternative='smaller'))
print()
# Гипотез о меньшем значении среднего перовой выборке
stat, p_value = ttest_ind(data1, data2, alternative = 'greater')
if p_value > alpha:
print('Принять гипотезу о меньшем значении среднего первой выборки')
else:
print('Отклонить гипотезу о меньшем значении среднего первой выборки')
print("Мощность критерия", smp.ttest_power(np.mean(data1)/np.std(data1),
nobs=len(data1), alpha = alpha, alternative='larger'))
for alpha in (0.1, 0.05, 0.01):
print('\nДоверительный интервал '+str(1 - alpha))
student_test(alpha, df.trestbps, df.thalach)
Доверительный интервал 0.9 Отклонить гипотезу о равенстве Мощность критерия 1.0 Отклонить гипотезу о большем значении среднего первой выборки Мощность критерия 0.0 Принять гипотезу о меньшем значении среднего первой выборки Мощность критерия 1.0 Доверительный интервал 0.95 Отклонить гипотезу о равенстве Мощность критерия 1.0 Отклонить гипотезу о большем значении среднего первой выборки Мощность критерия 0.0 Принять гипотезу о меньшем значении среднего первой выборки Мощность критерия 1.0 Доверительный интервал 0.99 Отклонить гипотезу о равенстве Мощность критерия 1.0 Отклонить гипотезу о большем значении среднего первой выборки Мощность критерия 0.0 Принять гипотезу о меньшем значении среднего первой выборки Мощность критерия 1.0
def w_test(alpha, data1, data2):
# Гипотеза о равенстве
stat, p_value = stats.mannwhitneyu(data1, data2, alternative='two-sided')
if p_value > alpha:
print('Принять гипотезу о равенстве')
else:
print('Отклонить гипотезу о равенстве')
# Гипотеза о большем значении среднего первой выборке'
stat, p_value = stats.mannwhitneyu(data1, data2, alternative='less')
if p_value > alpha:
print('Принять гипотезу о большем значении среднего первой выборки')
else:
print('Отклонить гипотезу о большем значении среднего первой выборки')
# Гипотез о меньшем значении среднего перовой выборке
stat, p_value = stats.mannwhitneyu(data1, data2, alternative='greater')
if p_value > alpha:
print('Принять гипотезу о меньшем значении среднего первой выборки')
else:
print('Отклонить гипотезу о меньшем значении среднего первой выборки')
for alpha in (0.1, 0.05, 0.01):
print('\nДоверительный интервал '+str(1 - alpha))
w_test(alpha, df.trestbps[:30], df.trestbps[-30:])
Доверительный интервал 0.9 Принять гипотезу о равенстве Принять гипотезу о большем значении среднего первой выборки Принять гипотезу о меньшем значении среднего первой выборки Доверительный интервал 0.95 Принять гипотезу о равенстве Принять гипотезу о большем значении среднего первой выборки Принять гипотезу о меньшем значении среднего первой выборки Доверительный интервал 0.99 Принять гипотезу о равенстве Принять гипотезу о большем значении среднего первой выборки Принять гипотезу о меньшем значении среднего первой выборки
sns.displot(x='trestbps', kde=True, data=df)
<seaborn.axisgrid.FacetGrid at 0x7fb474f45700>
sns.displot(x='thalach', kde=True, data=df)
<seaborn.axisgrid.FacetGrid at 0x7fb474e6b880>
def tests(alpha, data1, data2):
f = np.var(data1, ddof=1)/np.var(data2, ddof=1)
x = np.array(data1)
y = np.array(data2)
nun = x.size-1
dun = y.size-1
p_value = 1-stats.f.cdf(f, nun, dun)
print('\nТест Фишера:', end=' ')
if p_value > alpha:
print('Одно распределение')
else:
print('Разное рапределение')
stat, p_value = stats.levene(data1, data2, center='mean')
print('\nТест Левене:', end=' ')
if p_value > alpha:
print('Одно распределение')
else:
print('Разное рапределение')
stat, p_value = stats.bartlett(data1, data2)
print('\nТест Бартлета:', end=' ')
if p_value > alpha:
print('Одно распределение')
else:
print('Разное рапределение')
stat, p_value = stats.fligner(data1, data2, center='mean')
print('\nТест Флигнера-Килина:', end=' ')
if p_value > alpha:
print('Одно распределение')
else:
print('Разное рапределение')
for alpha in (0.1, 0.05, 0.01):
print('\nДоверительный интервал '+str(1 - alpha))
tests(alpha, df.trestbps[:30], df.thalach[:30])
Доверительный интервал 0.9 Тест Фишера: Одно распределение Тест Левене: Одно распределение Тест Бартлета: Одно распределение Тест Флигнера-Килина: Одно распределение Доверительный интервал 0.95 Тест Фишера: Одно распределение Тест Левене: Одно распределение Тест Бартлета: Одно распределение Тест Флигнера-Килина: Одно распределение Доверительный интервал 0.99 Тест Фишера: Одно распределение Тест Левене: Одно распределение Тест Бартлета: Одно распределение Тест Флигнера-Килина: Одно распределение
for c in ['pearson', 'spearman', 'kendall']:
plt.figure(figsize = (24, 12))
sns.heatmap(df.corr(method=c),
annot=True,
fmt=".2f")
plt.title(c)
plt.show()
con = pd.crosstab(df.sex, df.target)
ret = chi2_contingency(con)[1]
print(round(ret, 5))
0.0
ret = stats.fisher_exact(con)[1]
print(round(ret, 5))
0.0
print(mcnemar(np.array(con), exact=False)) #используем распределение хи квадрат
pvalue 0.00264477649712634 statistic 9.03763440860215
data = df[['sex', 'exang', 'target']]
CMH(data, 'sex', 'target', stratifier = 'exang')
df.corr()
| age | sex | cp | trestbps | chol | fbs | restecg | thalach | exang | oldpeak | slope | ca | thal | target | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| age | 1.000000 | -0.098447 | -0.068653 | 0.279351 | 0.213678 | 0.121308 | -0.116211 | -0.398522 | 0.096801 | 0.210013 | -0.168814 | 0.276326 | 0.068001 | -0.225439 |
| sex | -0.098447 | 1.000000 | -0.049353 | -0.056769 | -0.197912 | 0.045032 | -0.058196 | -0.044020 | 0.141664 | 0.096093 | -0.030711 | 0.118261 | 0.210041 | -0.280937 |
| cp | -0.068653 | -0.049353 | 1.000000 | 0.047608 | -0.076904 | 0.094444 | 0.044421 | 0.295762 | -0.394280 | -0.149230 | 0.119717 | -0.181053 | -0.161736 | 0.433798 |
| trestbps | 0.279351 | -0.056769 | 0.047608 | 1.000000 | 0.123174 | 0.177531 | -0.114103 | -0.046698 | 0.067616 | 0.193216 | -0.121475 | 0.101389 | 0.062210 | -0.144931 |
| chol | 0.213678 | -0.197912 | -0.076904 | 0.123174 | 1.000000 | 0.013294 | -0.151040 | -0.009940 | 0.067023 | 0.053952 | -0.004038 | 0.070511 | 0.098803 | -0.085239 |
| fbs | 0.121308 | 0.045032 | 0.094444 | 0.177531 | 0.013294 | 1.000000 | -0.084189 | -0.008567 | 0.025665 | 0.005747 | -0.059894 | 0.137979 | -0.032019 | -0.028046 |
| restecg | -0.116211 | -0.058196 | 0.044421 | -0.114103 | -0.151040 | -0.084189 | 1.000000 | 0.044123 | -0.070733 | -0.058770 | 0.093045 | -0.072042 | -0.011981 | 0.137230 |
| thalach | -0.398522 | -0.044020 | 0.295762 | -0.046698 | -0.009940 | -0.008567 | 0.044123 | 1.000000 | -0.378812 | -0.344187 | 0.386784 | -0.213177 | -0.096439 | 0.421741 |
| exang | 0.096801 | 0.141664 | -0.394280 | 0.067616 | 0.067023 | 0.025665 | -0.070733 | -0.378812 | 1.000000 | 0.288223 | -0.257748 | 0.115739 | 0.206754 | -0.436757 |
| oldpeak | 0.210013 | 0.096093 | -0.149230 | 0.193216 | 0.053952 | 0.005747 | -0.058770 | -0.344187 | 0.288223 | 1.000000 | -0.577537 | 0.222682 | 0.210244 | -0.430696 |
| slope | -0.168814 | -0.030711 | 0.119717 | -0.121475 | -0.004038 | -0.059894 | 0.093045 | 0.386784 | -0.257748 | -0.577537 | 1.000000 | -0.080155 | -0.104764 | 0.345877 |
| ca | 0.276326 | 0.118261 | -0.181053 | 0.101389 | 0.070511 | 0.137979 | -0.072042 | -0.213177 | 0.115739 | 0.222682 | -0.080155 | 1.000000 | 0.151832 | -0.391724 |
| thal | 0.068001 | 0.210041 | -0.161736 | 0.062210 | 0.098803 | -0.032019 | -0.011981 | -0.096439 | 0.206754 | 0.210244 | -0.104764 | 0.151832 | 1.000000 | -0.344029 |
| target | -0.225439 | -0.280937 | 0.433798 | -0.144931 | -0.085239 | -0.028046 | 0.137230 | 0.421741 | -0.436757 | -0.430696 | 0.345877 | -0.391724 | -0.344029 | 1.000000 |
y, x = dmatrices('target ~ age + sex + cp + trestbps + chol + fbs + restecg + thalach + exang + oldpeak + slope + ca + thal', data=df, return_type='dataframe')
vif = pd.DataFrame()
vif['VIF'] = [variance_inflation_factor(x.values, i) for i in range(x.shape[1])]
vif['variable'] = x.columns
vif
| VIF | variable | |
|---|---|---|
| 0 | 207.256646 | Intercept |
| 1 | 1.443474 | age |
| 2 | 1.161866 | sex |
| 3 | 1.284456 | cp |
| 4 | 1.170591 | trestbps |
| 5 | 1.150174 | chol |
| 6 | 1.087379 | fbs |
| 7 | 1.060998 | restecg |
| 8 | 1.613726 | thalach |
| 9 | 1.402001 | exang |
| 10 | 1.705857 | oldpeak |
| 11 | 1.642595 | slope |
| 12 | 1.202570 | ca |
| 13 | 1.147279 | thal |
for pr in ['age', 'sex', 'cp', 'trestbps', 'chol', 'fbs', 'restecg', 'thalach', 'exang', 'oldpeak', 'slope', 'ca', 'thal']:
aov = stats.f_oneway(df[pr], df.target)[1]
print(pr+' ~ target p-value =', aov)
age ~ target p-value = 0.0 sex ~ target p-value = 0.0004397191734875478 cp ~ target p-value = 2.845231775397828e-10 trestbps ~ target p-value = 0.0 chol ~ target p-value = 0.0 fbs ~ target p-value = 8.910856675010132e-27 restecg ~ target p-value = 0.6920275427852707 thalach ~ target p-value = 0.0 exang ~ target p-value = 4.686916079913056e-08 oldpeak ~ target p-value = 2.2270855709918155e-11 slope ~ target p-value = 2.970631656050962e-62 ca ~ target p-value = 0.004846478356075194 thal ~ target p-value = 4.408516795472436e-167
def regr_model(X, y):
X = np.array(X)[:, np.newaxis]
regr = LinearRegression()
cub = PolynomialFeatures(degree=3)
X_fit = np.arange(X.min(), X.max(), 1)[:, np.newaxis]
X_cub = cub.fit_transform(X)
regr = regr.fit(X, y)
y_lin_fit = regr.predict(X_fit)
regr = regr.fit(X_cub, y)
y_cub_fit = regr.predict(cub.fit_transform(X_fit))
plt.scatter(X, y)
plt.plot(X_fit, y_lin_fit,
label='linear',
lw=2)
plt.plot(X_fit, y_cub_fit,
label='cubic',
lw=2)
plt.legend()
plt.rcParams["figure.figsize"] = (6, 4)
regr_model(df.age, df.trestbps)